Self-organized chaos through polyhomeostatic optimization 
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The goal of polyhomeostatic control is to achieve a certain target distribution of behaviors, in 
contrast to homeostatic regulation which aims at stabilizing a steady-state dynamical state. We con- 
sider poly homeostasis for individual and networks of firing-rate neurons, adapting to achieve target 
distributions of firing rates maximizing information entropy. We show that any finite polyhomeo- 
static adaption rate destroys all attractors in Hopfield-like network setups, leading to intermittently 
bursting behavior and self-organized chaos. The importance of polyhomeostasis to adapting behav- 
ior in general is discussed. 



Introduction.- Homeostatic regulation plays a central 
role in all living, as well as in many technical applica- 
tions. Biological parameters, like the blood sugar level, 
the heart beating frequency or the average firing rates of 
neurons need to be maintained within certain ranges in 
order to guarantee survival. The same holds in the tech- 
nical regime for the rotation speed of engines and the 
velocity of airplanes, to give a few examples. 

Homeostasic control in the brain goes beyond the reg- 
ulation of scalar variables like the concentration of pro- 
teins and ions, involving the functional stability of neu- 
ral activity both on the individual as well as on a net- 
work level [D-[3|- We use here the term 'polyhomeosta- 
sis' for self-regulating processes aimed at stabilizing a 
certain target distribution of dynamical behaviors. Poly- 
homeostasis is an important concept used hitherto mostly 
implicitly and not yet well studied from the viewpoint 
of dynamical system theory. Polyhomeostasis is present 
whenever the goal of the autonomous control is the stabi- 
lization of a non-trivial distribution of dynamical states, 
polyhomeostatic control hence generalizes the concept of 
homeostasis. The behavior of animals on intermediate 
time scales, to give an example, may be regarded as 
polyhomeostatic, aiming at optimizing a distribution of 
qualitatively different rewards, like food, water and pro- 
tection; animals are not just trying to maximize a single 
scalar reward quantity. A concept loosely related to poly- 
homeostasis is homeokinesis, proposed in the context of 
closed- loop motion learning [4j, having the aim to sta- 
bilize non-trivial but steady-state movements of animals 
and robots. 

Here we study generic properties of dynamical systems 
governed by polyhomeostatic self-regulation using a pre- 
viously proposed model 5, 6] for regulating the firing-rate 
distribution of individual neurons based on information- 
theoretical principles. We show that polyhomeostatic 
regulation, aiming at stabilizing a specific target distribu- 
tion of neural activities gives rise to non-trivial dynami- 
cal states when recurrent interactions are introduced. We 
find, in particular, that the introduction of polyhomeo- 
static control to attractor networks leads to a destruction 
of all attractors resulting for large networks, as a function 
of the average firing rate, in either intermittent bursting 



behavior or self-organized chaos, with both states being 
globally attracting in their respective phase spaces. 

Firing-rate distributions - We consider a discrete-time, 
rate encoding artificial neuron with input x G [—00,00], 
output y € [0, 1] and a transfer function g(z), 



y(t+l) = g(a(t)x{t)+b(t)) , g(z) = 
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(1) 



The gain a(t) and the threshold —b(t)/a(t) in (JTJ arc 
slow variables, their time evolution being determined by 
polyhomeostatic considerations. 

Information is encoded in the brain through the firing 
states of neurons and it is therefore plausible to postu- 
late [5|, that polyhomeostatic adaption for the internal 
parameters a(t) and b(t) leads to a distribution p(y) for 
the firing rate striving to encode as much information 
as possible given the functional form (JTJ of the transfer 
function g{z). The normalized exponential distribution 
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maximizes the Shannon entropy [7j, viz the information 
content, on the interval y £ [0, 1], for a given expecta- 
tion value fi. A measure for the closeness of the two 
probability distributions p(y) and p\(y) is given by the 
Kullback-Leibler divergence [7[ 
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(3) 



which is, through (JTJ) , a function of the internal parame- 
ters a and b. The Kullback-Leibler divergence is strictly 
positive and vanishes only when the two distributions are 
identical. By minimizing D\(a,b) with respect to a and 
b one obtains Q the stochastic gradient rules 
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which have been called 'intrinsic plasticities' The re- 
spective learning rates e a and ej, are assumed to be small, 
viz the time evolution of the internal parameters a and b 
is slow compared to the evolution of both x and y. For 
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FIG. 1: (color online) Left: As a function of the average 
firing rate y,, the region of stability (shaded area) of the fix- 
point (y*,b*), see Eq. (0 of the one-neuron network in the 
quasistatic limit. Right: A self-coupling neuron with internal 
parameters a(t) and b(t). 



any externally given time series x(t) for the input, the 
adaption rules (UJ) will lead to a distribution of the out- 
put firing rates y(t) as close as possible, given the spec- 
ification ([T|) for the transfer function, to an exponential 
with given mean [i. 

Attractor relics.- fn deriving the stochastic gradient 
rules ([31) it has been assumed that the input x(t) is sta- 
tistically independent of the output y(t) (but not vice 
versa). This is not any more the case when a set of poly- 
homeostatically adapting neurons are mutually intercon- 
nected, forming a recurrent network, Here we will show 
that networks based on polyhomeostatic principles will 
generically show spontaneous and continuously ongoing 
activity. 

In a first step we analyze systematically the smallest 
possible network, viz the single-site loop, obtained by 
feeding the output back to the input, see Fig. [TJ a setup 
which has been termed in a different context 'autapse' 
Q. We use the balanced substitution x — > y — 1/2 in 
Eqs. ([1} and ((3]), the complete set of evolution rules for 
the dynamical variables y(t), a(t) and bit) is then 

y(t + l) = g(a(t)[y(t)- 1/2} +b(t)) 
b(t+l) = b(t) + e b Ab(t) (5) 
a(t + 1) = a(t) + e a (l/a(t) + [y{t) - l/2]A6(t)' 



with 



Ab(t) = l-(2 + \)y(t + l) + \y 2 {t 



(6) 



Note, that Ab(t) in ([6]) depends on y(t + 1), and not 
on y(t), as one can easily verify when going through the 
derivation of the rules ^ for the intrinsic plasticity. The 
evolution equations (|5|) are invariant under y <-> (1 - y), 
b o (—5), a -H> a and A «-» (—A), the later corresponding 
to the interchange of n «-> (1 — fi). 

We first consider the quasistatic limit e a <C viz 
a(t) ~ a is approximatively constant. The fixpoint 
(y*, b*) in the (y, b) plane is then determined by 
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FIG. 2: (color online) The time dependence of y(t) (red stars) 
b(t) (solid line) and of a(t) — 4 (dashed line) for the balanced 
one-site problem e a = = e = 0.01 and fi = 0.28; the 
horizontal dotted lines are guides to the eye. The gain a(t) is 
initially small and the system relaxes, since e < 1, fast to a 
fixpoint of y = g(a[y — 1/2] +6). Once a(t) surpasses a certain 
threshold, compare Fig. [TJ the fixpoint becomes unstable and 
the system starts to spike spontaneously. 



A straightforward linear stability analysis shows, that 
the fixpoint (y* , b*) remains stable for small gains a and 
becomes unstable for large gains, see Fig. [T] We now 
go beyond the quasistatic approximation and consider 
a small but finite polyhomeostatic adaption rate e a for 
the gain. Starting with a small gain a we see, compare 
Eq. <j5j) , that the gain necessarily grows until it hits the 
boundary towards instability; for small Ab the growths 
of the gain is a(t) ~ \/t. 

In other words, a finite adaption rate e a for the gain 
turns the fixpoint attractor (y*,b*) into an attractor 
relic and the resulting dynamics becomes non trivial and 
autonomously self-sustained. This route towards au- 
tonomous neural dynamics is actually an instance of a 
more general principle. Starting with an attractor dy- 
namical system, viz with a system governed by a finite 
number of stable fixpoint, one can quite generally turn 
these fixpoints into unstable attractor ruins by coupling 
locally to slow degree of freedoms [Tjl [ll| , in our case the 
slow local variables are a(t) and b(t). Interestingly, there 
are no saddlepoints present in attractor relic networks in 
general, and in ([5]) in particular, and hence no unstable 
heteroclines, as in a proposed alternative route to tran- 
sient state dynamics via heteroclinic switching. 12|, [13 [ 

In Fig.[5]the time evolution is illustrated for A = 3.017, 
which corresponds to // = 0.28, see Eq. and e a = 
£5 = e = 0.01. The system remains in the quasistation- 
ary initial regime until the gain a surpasses a certain 
threshold. The initial quasistationary fixpoint becomes 



therefore unstable via the same mechanism discussed an- 
alytically above, see Eq. (J7J, for the regime e Q — > 0, com- 
pare Fig. [T] The output activity y(t) oscillates fast be- 
tween two transient fixpoints of y = g{a[y — 1/2] + b), 
having a high and a low value respectively. This spiking 
behavior of the neural activity is driven by spontaneous 
oscillations in the threshold —b(t)/a(t), shifting the in- 
tersection of g(a[y — 1/2] + b) with y forth and back. 



3 



+ 


2- 








1 - 








0- 






>> 


■4 



a(t)-4 

. . . 

yd) 



b(t)+2 



1000 



2000 



3000 



FIG. 3: (color online) The time dependence of y(t) (red stars) 
b(t) + 2 (solid line) and of a(t) — 4 (dashed line) for the one- 
site problem with inhibitory self-coupling x — s> 1/2 — y, e a = 
et = e = 0.01 and = 0.28. A Hopf-bifurcation occurs for the 
output y(t) when the initial quasistationary fixpoint becomes 
unstable, giving place to a new fixpoint of period two. 



Evaluating the local Lyapunov exponents we find that 
the trajectory is stable against perturbations for the tran- 
sient states close to one of the transient fixpoints of 
y = g(a[y — 1/2] + b) and sensitive to perturbations and 
external modulation during the fast transition periods, 
an instantiation of the general notion of transient state 
dynamics El • 

Inhibitory self- coupling - So far we discussed, see §5§ 
and Fig.rjJ a neuron having its output y(t) coupled back 
excitatorily to its input via x — > y — 1/2. The dynam- 
ics changes qualitatively for an inhibitory self-coupling 
x — > 1/2 — y, see Fig. [3J There is now only a single 
intersection of g(a[—y + 1/2] + b) with y. This intersec- 
tion corresponds to a stable fixpoint for small gains a. A 
Hopf-bifurcation [3] occurs when a(t) exceeds a certain 
threshold and a new fixpoint of period two becomes sta- 
ble. The coordinates of this fixpoint of period two slowly 
drift, compare Fig. [3J due to the residual changes in a(t) 
and b(t). Interestingly, the dynamics remains non-trivial, 
as a consequence of the continuous adaption, even in the 
case of inhibitory self-coupling. 

Self organized chaos- We have studied numerically 
fully connected networks of i — 1,..,N polyhomeostat- 
ically adapting neurons (@|, coupled via 



(8) 



The synaptic strengths are Wij = ±l/\/N — 1, with in- 
hibition/excitation drawn randomly with equal proba- 
bility. The adaption rates are e a = = 0.01. We con- 
sider homogeneous networks where all neurons have the 
identical fj, for the target output distributions ([2]), with 
fj, = 0.15, 0.28 and fi = 0.5. 

The activity patterns presented in the inset of Fig. @] 
are chaotic for the /i = 0.28 network and laminar with in- 
termittent bursting for the network with /i = 0.15. This 
behavior is typical for a wide range of network geome- 
tries, distributions of the synaptic weights and for all ini- 
tial conditions sampled. The respective Kullback-Leibler 
divergences D\ decrease with time passing, as shown in 
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FIG. 4: (color online) Mean Kullback-Leibler divergence D\, 
for N = 500 networks, of the real and target output dis- 
tributions as a function of number of iterations t. The solid 
black, dashed-dotted red and the dashed blue lines correspond 
respectively to target firing rates fi = 0.15, /j, = 0.28, and 
fi = 0.5. Inset: Output activity of two randomly chosen neu- 
rons for target average firing rates n = 0.28 (top) and n — 0.15 
(bottom). 



Fig. 01 due to the ongoing adaption process (jj). D\ be- 
comes very small, though still finite, for long times in the 
self-organized chaotic regime (/x = 0.28, 0.5), remaining 
substantially larger in the intermittent bursting phase 
(fx = 0.15). 

We have evaluated the global Lyapunov exponent, [l5| 
finding that two initial close orbits diverge until their 
distance in phase space corresponds, within the available 
phase space, to the distance of two randomly selected 
states, the tellmark sign of deterministic chaos. [3] The 
corresponding global Lyapunov exponent is about 5% 
per time-step for fi = 0.5 and initially increases with 
the decrease of [i, until the intermittent-bursting regime 
emerges, after which the global Lyapunov exponent de- 
clines with the decrease of fx. 

The system enters the chaotic regime, in close anal- 
ogy to the one-side problem discussed previously, when- 
ever the adaptive dynamics (|4|) has pushed the individ- 
ual gains ai(t), of a sufficient number of neurons, above 
their respective critical values. Hence, the chaotic state 
is self-organized. Chaotic dynamics is also observed in 
the non-adapting limit, with e a ,£b — ► 0, whenever the 
static values of are above the critical value, in agree- 
ment with the results of a large- TV mean field anal ysis 
of an analogous continuous time Hopfield network. [16] 
Subcritical static dj lead on the other side to regular dy- 
namics controlled by point attractors. 

In Fig. [5] we present the distribution of the output ac- 
tivities for networks with target mean firing rates [i — 
0.28 and fj, — 0.5. Shown are in both cases the firing rate 
distributions p{y) of the two neurons having respectively 
the largest and the smallest Kullback-Leibler divergence 
([3j) with respect to the target exponential distributions 
@. Also shown in Fig. [S] are the distributions of the 
network-averaged output activities. 
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FIG. 5: (color online) Output distributions of the two neu- 
rons with highest (blue diamonds) and lowest (red circles) 
Kullback-Leibler divergence © compared to the mean out- 
put distribution (dashed green line) and the target exponen- 
tial output distribution (full black line). For a network with 
N = 500 neurons and a target output mean firing rate (A) 
fi — 0.28 and (B) li = 0.5 (identical for all neurons). Insets: 
For the same parameters the mean output distribution of the 
single self-coupled neuron, see Fig. [1] 



mittency to chaos as a function of /x and system sizes will 
be given elsewhere. 

Conclusions. -It is well known 



17], that individual and 



networks of spiking neurons may show bursting behav- 
ior. Here we showed, that networks of rate encoding 
neurons are intermittently bursting for low average fir- 
ing rates, when polyhomeostatically adapting, entering a 
fully chaotic regime for larger average activity level. Au- 
tonomous and polyhomeostatic self-control of dynamical 
systems may hence lead quite in general to non-trivial 
dynamical states and novel phenomena. 

Polyhomeostatic optimization is of relevance in a range 
of fields. Here we have discussed it in the framework of 
dynamical system theory. Alternative examples are re- 
source allocation problems, for which a given limited re- 
source, like time, needs to be allocated to a set of uses, 
with the target distribution being the percentages of re- 
source allocated to the respective uses. These ramifica- 
tions of polyhomeostatic optimization will be discussed 
in further studies. 
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The data presented in Fig.[5]shows, that the polyhome- 
ostatically self-organized state of the network results in 
firing-rate distributions close to the target distribution. 
This result is quite remarkable. The polyhomeostatic 
adaption rules (U) are local, viz every neuron adapts its 
internal parameters aj and hi independently on its own. 
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both regimes. A qualitative difference can be observed 
when considering the subspace of activities {y\, . . . , Un)- 
Evaluating the global Lyapunov exponent in this sub- 
space, viz the relative time evolution (Ayi, . . . , Ay^) of 
differences in activities, we find this restricted global Lya- 
punov to be negative in the intermittent and positive in 
the chaotic regime. Details on the evolution from inter- 
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